A technique for improving petroleum products forecasts using grey convolution models and genetic algorithms

Forecasting energy consumption is a major concern for policymakers, oil industry companies, and many other associated businesses. Though there exist many forecasting tool, selecting the most appropriate one is critical. GM(1,1) has proven to be one of the most successful forecasting tool. GM(1,1) does not require any specific information and can be adapted to predict energy consumption using a minimum of four observations. Unfortunately, GM(1,1) on its own will generate too large forecast errors because it performs well only when data follow an exponential trend and should be implemented in a political-socio-economic free environment. To reduce these short-comings, this paper proposes a new GM(1,n) convolution model optimized by genetic algorithms integrating a sequential selection mechanism and arc consistency, abbreviated Sequential-GMC(1,n)-GA. The new model, like some recent hybrid versions, is robust and reliable, with MAPE of 1.44%, and RMSE of 0.833.• Modification, extension and optimization of grey multivariate model is done.• The model is very generic can be applied to a wide variety of energy sectors.• The new hybrid model is a valid forecasting tool that can be used to track the growth of households’ energy demand.


Introduction
In 2017, the household sector was Cameroon's largest energy consumer, accounting for 70% of overall energy consumption, thereby over distancing transports (15%) and industries (11%) [1] . It was also the fastest growing sector over the period 2001-2010, going from 1942.78 kilotons of oil equivalent (ktoe) to 3995.57 ktoe, although strongly dominated by biomass with a share of 94.74% (3785.6 ktoe), followed by petroleum products (PP) (125.13 ktoe, 3.13%) and electricity (84.84 ktoe, 2.12%) [1] . In the following years, this trend is projected to continue, due to major urbanization projects initiated by the State of Cameroon. Cameroon thus finds itself torn between two major challenges: ensuring the development of its household sector without ignoring climatic issues. To achieve these projects, the certainty on energy consumption prediction is quite important.

Interest of study
The legal obligation imposed by regulatory authorities is one reason for conducting forecasting studies in Cameroon's oil market. Besides, in cases where distributors cause an imbalance in the supply chain of PPs, related regulations and laws impose a high penalty payment on them [2] . This has the effect of putting pressure on distributors to accurately forecast end-users' consumption.
Forecast studies are also important because PPs supply is provided through spot markets, where the bulk of the demand is in urban areas. Most households in Cameroon consume the bare minimum of their PP needs, as the latter are only used for lighting (kerosene) and/or cooking (kerosene and liquefied petroleum gas (LPG)) [2] . It is therefore incumbent on the State, and according to the market's structure, to manage its own strategic stocks in order to regulate PP's consumption.

Related works
There exist many forecasting techniques [3] , however, genetic algorithms (GA) and grey models (GM) are different from all others by their simplicity. In addition, they can be applied with a minimum of four data points. Unfortunately, GMs alone are not sufficient to forecast energy demand because, very often, forecast results are not sufficiently precise and sometimes, too strong assumptions are made which could irremediably lead to unacceptable forecast errors [4] . GA have gained the attention, as they are robust stochastic search algorithms that are used to solve various problems, including forecasting. GA have several considerable advantages over conventional methods. They can provide practical solutions, by analysing the search space from many starting points, without making any prior assumption about the model or the underlying function [5] .
Amongst existing grey models, GM(1,1) has had success in many fields including energy, earthquakes, finance, food production, education, road safety and transports. However, GM(1,1) works well when data follow an exponential trend and this is only applicable in a political-socio-economic free environment. Yet, such a neutral environment does not reflect the reality of a generalized grey system. This is why data collected from such an environment cannot be modelled with GM(1,1) because there would be incomplete information. To complete this information, more variables must be inserted into the system. This tends to suggest that a multivariate grey model (GM(1,n)) is a way forward to whiten the system. Unfortunately, as demonstrated by Tien [5] , GM(1,n) and GM(m,n) all lead to inaccurate modelling because of their structural defects.
The traditional GM(1,n) can be improved by adjusting its structure to make it more flexible either by directly optimizing and stabilizing its parameters, and by considering the hysteresis of the related variables [6] . Other lines of thought include new accumulation generating operation, hybridization, or by integrating a functional mechanism into the GM's structure. To date, no work has been carried out on the integration of specific mechanisms into the forecasting framework of GM(1,n). This could take into consideration key drivers of energy consumption, by connecting GM with GA.
This method proposes a sequential grey multivariable convolution model optimized by GA (Sequential-GMC(1,n)-GA). The technical innovation here is to integrate the sequential selection mechanisms in GMC(1,n) and the arc-consistency (AC) in GA as a filter with the aim of reducing CPU execution time. In this perspective, GMC(1,n) is trained by GA in order to enhance all parameters and improve prediction accuracy. The motivation of combining GMC(1,n) with GA are threefold, namely: To whiten the grey system under consideration rendering its information completely extractable. When this is done, the hybrid model can then correctly capture predictable structure in the demand history including: trend; seasonality; and special events (e.g. peaks) that could impact demand. To optimize all GM parameters rather than using theoretical values. This ensures that all possible solutions have been explored and tested, given way to unbiased predictions, meaning that the estimated forecast is not projecting too high or too low. To gather the most recent data characteristics which minimizes predicting errors, CPU execution time, and the amount of retracing required because inescapable faults may be quickly identified.
As novelty, this method proposes a GMC that fully excavates the evolution law of multivariate time series without modifying the model's framework. Sequential-GMC(1,n)-GA proactively improves all model parameters in each projected time-frame, regardless of data patterns. Thus, Sequential-GMC(1,n)-GA eliminates the inconsistency issue between parameterization approach that is based on grey derivative and the model's minimal sum of squares of prediction errors, resulting in more accurate forecasts. Finally, the model is capable of using input data with different sizes and still compete with performant forecasting models. Section 2 that follows presents the method. Section 3 presents the application using real data. Setup for model's evaluation is given in Section 4, while Section 5 is the conclusion.
Step 3. Designing the background value of the system The system's background value is defined as in Eq. (4) : The horizontal adjustment coefficient is taken such that 0 < < 1 [7] . The value of is generally taken as 0.5 but it should be chosen so as to reduce the forecast errors [7] .
Step 4. Establishing Grey systems of equations Grey system is trained by establishing a link between known and unknown sequences: is the development coefficient, is GMC(1,n) parameter, while =1 , …, are grey input coefficients. Tien [5] starts by considering the right hand side of Eq. (5) as a function ( ) . Then, using the trapezoid formula, and integrating both sides of Eq. (5) from − 1 to , Eq. (5) is approximated by the following difference equation: As a result, Eq. (6) can be regarded as a linear equation system in relation to the coefficients [ 2 3 … ] . These coefficients are calculated by least squares method. Thus, by applying least squares, ] are determined as shown in Eq. (7) : Step 5. Solving the system's differential equation The solution to Eq. (5) is [5] : Whose value can only be approximated with numerical methods because of the presence of convolution integral. With trapezoid formula, Eq. (8) becomes: Recall that ̂ (1) 1 (1) = (0) 1 (1) , and ℎ ( ) in Eq. (9) is defined as: ℎ ( ) = 0 if < 2 ; else ℎ ( ) = 1 Step 6. Inverse AGO (IAGO) Finally, forecasted values of ̂ (0) 1 ( ) are obtained by IAGO as shown in Eq. (10) : Optimizing GMC(1,n) with GA using arc consistency (AC) technique The parameters , , and =1 , …, affect the precision of GMs. In this method, the best values of these parameters are determined experimentally (instead of using arbitrary values reported in the literature).

Formalism of constraint satisfaction problems (CSPs)
A The arity of a constraint is the number of variables on which the constraint relates. A CSP with arity larger than 2 is an -ary CSP. A solution is the assignment of a value to all variables of the problem such that each constraint is satisfied. CSP is consistent if and only if it admits a solution. A tuple of the constraint  is said to be allowed if and only if ∈   and said to be supported for a value ∈  , ∈ (  ) , if and only if is allowed and contains in the position corresponding to in the constraint. Checking whether a given tuple is allowed by a constraint is called consistency testing. The possible combinations of variable instantiations in a CSP create a search space that can be seen as a search tree. We used AC algorithms to non-binary CSPs with nAC3 -ary algorithm based on the principle of AC3 [8] .
For constraint checking, we used the more popular nFC3 forward search algorithm. When instantiating a variable, nFC3 removes from the current domains of future variables all values incompatible with this instantiation. When a new variable is considered, we can then be sure that all the values of its current domain are consistent with the past variables. The computational steps of nAC3 -ary algorithm based on the principle of AC3 is given below: CSP Algorithm

Start
Step 1 : Step 3 : Create a constraint set with and  after considering the constraints 
where is the number in decimal form which is represented in binary form, is the number of bits, and are the lower and upper threshold values respectively. Selection, crossover and mutation operators are then used (over a series of generations) by a standard GA to guide ( ) towards convergence at the global optimum [9] . Once selection, reproduction and mutation operators have been applied, new phenotypes (population) are created (which theoretically has a better forecast accuracy). The generational counter is incremented by one. The whole process is reiterated until is reached or until the desired forecast accuracy is achieved.
Unluckily, GA are well-known for being CPU-intensive (they sometimes take a long time when executing a given problem). To tackle this issue, we assumed that each individual in ( ) is a CSP that we are attempting to solve using AC, with ( = 1 ) serving as the outcome. Further specifics are in [10] .

Sequential GMC(1,n) prediction model optimized by GA
, , and =1 , …, are found in the literature to be 0.5 and 4 for and respectively, while and =1 , …, are calculated only once during the modelling stage. However, periodic values of these parameters can significantly increase accuracy of GMC(1,n) [5] . Hence, instead of using values reported in the literature, we rather calculate them experimentally (using AC-GA) for each forecast period ℎ .
Sequential mechanism makes use of the most up-to-date simulations for prediction by depicting the latest characteristics of data. This allows the sequential mechanism to improve forecasting accuracy. This mechanism employs sequences to model GMC(1,n), and consumption data for forecasting. At each loop, new values of , , , =1 , …, are computed and optimized by GA. The process is described below:  This process is called sequential GMC(1,n) prediction optimized by GA (sequential-GMC(1,n)-GA) and is summarized in Fig. 1 . The above-mentioned instructions have been written with the AC (steps ii and iv). It is therefore a Sequential-GMC(1,n)-GA with AC. By editing the instruction "with AC" in steps ii and iv, it becomes a Sequential-GMC(1,n)-GA without AC.

Application of Sequential-GMC(1,n)-GA using real data
Useful variables must be separated from the unnecessary ones. Stepwise regression analysis (SRA) defines predictor variables that most accurately characterize the variable to be predicted [11] . For this, we start by regressing the dependant variable on a single independent variable and set the significance level and values of r (F-to-remove) and e (F-to-enter); being the Fisher statistic. The statistical significance of the F-test and the decrease in the sum of the squared error constitute the benchmark for adding or eliminating a variable. After inserting a new variable in the regression model, the partial F-value is calculated and compared to e and r ; if > e , then we include this variable; otherwise, if < r , it is eliminated and will never return to the model again. All useful variables used for this work are in annual frequency and are presented in Table 1 .
The modelling period was selected for 1992-2008, while the test period is 2009-2017. Dataset is divided into simulation (or training) set and validation (or test) set to ensure that the model is neither overfitting nor underfitting. Training dataset spans from 1992 to 2008 while validation dataset is from 2009 to 2017.
Prices and historical consumption were collected from HPSF ( https://www.csph.cm ) and confirmed by CCPD data ( https://www.scdp.cm ). Urbanization rate and number of households were provided by the National Institute of Statistics ( https://www.minefop.gov.cm ), while real income was collected from the World Bank statistics ( https://data.worldbank.org ).

Performance measurements
Threshold levels for MAPE, and 2 are given in Table 2 . For RMSE and AAE, low values close to 0 indicate the best predictions [12] . The speed of convergence and stability, i.e. quality and reliability, are used to assess GA's performance. The best solution obtained in relation to the number of cycles (computational cost) is used to determine speed of convergence, while the standard deviation of the results from 10 optimization runs per algorithm is used to determine stability.

Setup for model's evaluation
Sequential-GMC(1,n)-GA is initialized with parameters as indicated in Table 3 . Once approximate values of , , , =1 , …, are obtained, they are encoded into binary digits and integrated into GA. Amongst these parameters, mutation probability must be very low, as low as 0.05 or even smaller. Because a higher value could destroy the solution.
is set as the inverse of chromosome length. Crossover probability depends on the problem at hand. Its optimal value is estimated after different runs based on the targeted Table 2 Accuracy measures and threshold levels.
mean absolute percentage error, * coefficient of correlation, * * adjusted . † root mean squared error,. † † average absolute error.  precision . If optimal values of and are settled, then number of generations and population size are simply set arbitrarily. However, a population size of 100 and 50 generations is advised.

Simulations and model validation
The optimal values of , , , =1 , …, are presented in Table 4 and are different from theoretical values reported in the literature. We note that without AC, execution takes 97 s and consistency is achieved at the 34th generation. Thereafter, each chromosome in ( ) is assumed to be a CSP while each gene is considered as a variable. The global optimum is obtained at the 13th generation after generating the CSP at random and removing inconsistent values from the variable domains. The new execution takes 5 s. So, AC technique reduced execution time by 92 s. This gain also explains from 34 to 13 generations. As a result, GA needed fewer iterations to converge to a global optimum. AAEs for each model are shown in Figs. 2 (a-c) and 3 (a-c). Prediction deviations given by Sequential-GMC(1,n)-GA are insignificant compared to GMC(1,n) and OGMC(1,n). Graphs of real and predicted consumption are shown in Figs. 4 and 5 , with annual residual

Convergence and stability tests
We ran 10 trials of Sequential-GMC(1,n)-GA with and without AC, and the average value and standard deviation was computed to compare the speed of convergence and stability of the method (see Figs. 6-8 ).   Using LPG data, after 50 evaluations of Sequential-GMC(1,n), GA without AC finds mean solution at the 21st evaluation with MAPE of 4.622%. Using AC techniques combined with GA, mean solution is achieved at the 11th generation with MAPE of 1.178% ( Fig. 6 ). Sequential-GMC(1,n)-GA without AC exhibit early convergence, and Sequential-GMC(1,n)-GA with AC catches up at around 7 evaluations. This better converging algorithm also display high stability (see Fig. 8 ).
With kerosene data, after 50 evaluations, Sequential-GMC(1,n)-GA without AC finds a solution at the 30th evaluation with MAPE of 1.4% whereas Sequential-GMC(1,n)-GA with AC technique finds a solution at the 13th evaluation with 0.388% MAPE (see Fig. 7 ). As for stability, Fig. 8 reveals a large difference between Sequential-GMC(1,n)-GA with AC and Sequential-GMC(1,n)-GA without AC. It is already well known that GA slows down its convergence speed if it needs to fine-tune weights and biases. However, there are two possible explanations why Sequential-GMC(1,n)-GA with AC outperforms Sequential-GMC(1,n)-GA without AC. First, the threshold to stop GA without AC is too low for this task so that efficiency was lost when a plateau is encountered. The second reason is that Sequential-GMC(1,n)-GA with AC is faster even in the initial search of near optimal solutions.

Conclusion
In this paper, we proposed a Sequential-GMC(1,n)-GA hybrid model to forecast households' PP consumption. Forecasting models were produced with training datasets from 1992 to 2008, while validation datasets include data for the period 2009-2017. The results show that Sequential-GMC(1,n)-GA hybrid model is suit to predict energy demand than competing models. This superiority results from two facts: first, instead of using theoretical values of and , all parameters ( , , , =1 , …, ) are determined experimentally and optimized before each forecasting period. Second, the use of a sequential method and AC gathers the most recent data characteristics, reducing forecasting errors and CPU execution time. As a result, this study demonstrates that GM's performance can still be improved. Overall, the model has good adaptability and feasibility, and can extract a grey system's evolution law.
For forecasting purposes, the proposed model can be chosen out of other competing models provided the time series do not have a high degree of periodicity, volatility or are free from price shocks. This is because the Sequential-GMC(1,n)-GA will be unable to completely extract the evolution law due to the presence of disturbing characteristics. However, this issue can be solved by including a term that accounts for nonlinearities in a time-varying GMC(1,n) or by performing data preprocessing.

Funding
This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data Availability
Data will be made available on request.